---
title: 'Replication for "Thin-skinned Leaders: Regime Legitimation, Protest Issues
  and Repression in Autocracies"'
subtitle: 'Analysis script'  
author: "Eda Keremoglu, Sebastian Hellmeier, Nils B. Weidmann"
date: "2020"
output: html_document
---

```{r setup, include=FALSE, warning=FALSE}
knitr::opts_chunk$set(echo = T, include = T, warning = F)
```

```{r}
# load packages
library(tidyverse)
library(countrycode)
library(survival)
library(lubridate)
library(multiwayvcov)
library(stargazer)
library(lmtest)
library(sandwich)
library(xtable)

source("interaction_plot.R")

options(scipen = 99999, max.print = 99999)
```

load data
```{r}
load("data.Rda")
```

Figure 1

```{r, echo=FALSE, include=T, eval=T}

# select vdem claims data and arrange for plotting
load("vdem_claims.Rda")

vdem_claims <- vdem_claims %>% 
  mutate(country_name = countrycode(cowcode, "cown", "country.name"), year = as.numeric(as.character(year))) %>% 
  tidyr::pivot_longer(-c(country_name, year)) %>%
  filter((country_name == "Venezuela" | country_name == "Cuba") &
           year > 1999) %>% 
  dplyr::mutate(claims = stringr::str_remove(name, "_codelow|_codehigh"),
                    type = ifelse(stringr::str_detect(name, "_codelow"), "codelow", "value"),
                    type = ifelse(stringr::str_detect(name, "_codehigh"), "codehigh", type)) %>%
  tidyr::pivot_wider(names_from = (type),
                         values_from = value,
                         id_cols = c(country_name, year, claims)) %>%
  tidyr::drop_na()


# plot claims for Venezuela
ggplot(data = vdem_claims[vdem_claims$country_name=="Venezuela",],
       aes(x = year, y = value, group = claims, linetype = claims)) + 
  geom_line()  + 
  geom_ribbon(aes(x = year, ymin = codelow,
                  ymax = codehigh), alpha = .07) +
  xlab("") + ylab("Claim strength") + ylim(0,4) +
  scale_linetype_manual(name = "", values = c(2, 1, 3, 4),
                        labels = c("Ideology", "Leader", "Performance", "Rational-legal")) +
  scale_x_continuous(breaks = seq(2000, 2020, 2)) +
  theme(panel.background = element_rect(fill = "white", colour = "black", size = 0.5, linetype = "solid"),
        axis.text.x = element_text(angle = 45, hjust = 1, size = 10),
        panel.grid.major.x = element_line( size=.05, color="grey"),
        legend.position = "bottom") 


# plot claims for Cuba
ggplot(data = vdem_claims[vdem_claims$country_name=="Cuba",],
       aes(x = year, y = value, group = claims, linetype = claims)) + 
  geom_line()  + 
  geom_ribbon(aes(x = year, ymin = codelow,
                  ymax = codehigh), alpha = .07) +
  xlab("") + ylab("Claim strength") + ylim(0,4) +
  scale_linetype_manual(name = "", values = c(2, 1, 3, 4),
                        labels = c("Ideology", "Leader", "Performance", "Rational-legal")) +
  scale_x_continuous(breaks = seq(2000, 2020, 2)) +
  theme(panel.background = element_rect(fill = "white", colour = "black", size = 0.5, linetype = "solid"),
        axis.text.x = element_text(angle = 45, hjust = 1, size = 10),
        panel.grid.major.x = element_line( size=.05, color="grey"),
        legend.position = "bottom") 

```

Table 1

```{r, echo=FALSE, include=T, eval=T}
# load summary from BTM model
summary_LDA <- read_csv("summary_BTM.csv")

summary_LDA
```

Figure 2

```{r, echo=FALSE, include=T, eval=T}
# load btm model
load("model.Rds")

# take topic frequency from BTM model
topics_thetas <- model_BTM$theta

# add variable for topic labels 
topics <-  c("Background","Repression","Governance",
             "Incumbents/Officials",
             "Elections","Opposition/Repression","Regime",
             "Economy/Living conditions")

# combine both
topic_freqs <-cbind.data.frame(topics, topics_thetas)

  # plot topic intensity as bar plot
ggplot(topic_freqs[topic_freqs$topics!="Background",], aes(x=reorder(topics, -topics_thetas),
                                                           y=topics_thetas)) +
  geom_bar(stat = "identity", color = "black", fill="lightgray",width = 0.5) + ylab("Topic frequency") +
  xlab("") + 
  theme(panel.background = element_rect(fill = "white", colour = "black", size = 0.5, linetype = "solid"),
        axis.text.x = element_text(angle = 45, hjust = 1, size = 10))
```

Figure 3a and 3b

```{r, echo=FALSE, include=T, eval=T}
# Venezuela 2014
venezuela_2014 <- data %>% 
  ungroup() %>% 
  mutate(year = str_sub(event_date, 1, 4) , 
         week = floor_date(event_date, unit = "weeks",
                           week_start = getOption("lubridate.week.start", 1))) %>% 
  filter(cowcode==101 & year==2014) %>% 
  select(event_id, event_date, week, topic2:topic8) %>% 
  group_by(week) %>% 
  summarise(elections = sum(topic5, na.rm = T) / length(unique(event_id)),
            economy = sum(topic8, na.rm = T) / length(unique(event_id)),
            regime_incumb = (sum(topic7, na.rm = T) + sum(topic4, na.rm = T)) / length(unique(event_id))) %>% 
  mutate(other = 1 - (elections + economy + regime_incumb)) %>% 
  gather(., "economy", "regime_incumb", "elections", "other",
         key = "topic", value = "intensity") %>% 
  ungroup() %>% 
  mutate(topic = factor(topic, levels=c("elections","economy", "regime_incumb", "other")))

# Armenia 2013
armenia_2013 <- data %>% 
  ungroup() %>% 
  mutate(year = str_sub(event_date, 1, 4),
         week = floor_date(event_date, unit = "weeks",
                           week_start = getOption("lubridate.week.start", 1))) %>% 
  filter(cowcode==371 & year==2013) %>% 
  select(event_id, event_date, week, topic2:topic8) %>% 
  group_by(week) %>% 
  summarise(elections = sum(topic5, na.rm = T) / length(unique(event_id)),
            economy = sum(topic8, na.rm = T) / length(unique(event_id)),
            regime_incumb = (sum(topic7, na.rm = T) + sum(topic4, na.rm = T)) / length(unique(event_id))) %>% 
  mutate(other = 1 - (elections + economy + regime_incumb)) %>% 
  gather(., "economy", "regime_incumb", "elections", "other",
         key = "topic", value = "intensity") %>% 
  ungroup() %>% 
  mutate(topic = factor(topic, levels=c("elections","economy", "regime_incumb", "other")))

# plot Venezuela
ggplot(data = venezuela_2014, aes(x=week, y=intensity,  fill=topic)) + 
  geom_bar(stat = "identity", position = position_stack(reverse = TRUE), width=6) +
  theme_classic() +
  xlab("") + ylab("Topic probability") +
  scale_fill_manual(name = "", labels = c("Elections", "Economy", "Regime", "Other"),
                  values = c("#525252", "#cccccc", "#969696",  "#f7f7f7")) +
  theme(legend.position="bottom",
        legend.text=element_text(size = 11),
        axis.text=element_text(size = 11),
        axis.title.y = element_text(size = 12)) +
  scale_x_date(limits = as.Date(c('2014-01-01','2014-12-31')))

# plot Armenia
ggplot(data = armenia_2013, aes(x=week, y=intensity,  fill=topic)) + 
  geom_bar(stat = "identity", position = position_stack(reverse = TRUE), width = 6) +
  theme_classic() +
  xlab("") + ylab("Topic probability") +
  scale_fill_manual(name = "", labels = c("Elections", "Economy", "Regime", "Other"),
                  values = c("#525252", "#cccccc", "#969696",  "#f7f7f7")) +
  theme(legend.position="bottom",
        legend.text=element_text(size = 11),
        axis.text=element_text(size = 11),
        axis.title.y = element_text(size = 12)) +
  scale_x_date(limits = as.Date(c('2013-01-01','2013-12-31')))
```

Table 2

```{r}
# issues with fe
m1 <- clogit(repression_binary ~ issue_performance + issue_rat_legal + issue_leader + issue_regime + topic2 + 
               year + strata(cowcode) + cluster(cowcode), data = data, method = "efron")

# with claims and fe
m2 <- clogit(repression_binary ~ issue_performance + issue_rat_legal + issue_leader + issue_regime + topic2 +
            v2exl_legitperf + v2exl_legitratio  + v2exl_legitlead + 
            year + strata(cowcode) + cluster(cowcode), data = data, method = "efron")

# with interactions
m3 <- clogit(repression_binary ~ issue_performance*v2exl_legitperf + issue_rat_legal*v2exl_legitratio + issue_leader*v2exl_legitlead + issue_regime + topic2 +
            year + strata(cowcode) + cluster(cowcode), data = data, method = "efron")

# with event-level controls
m4 <- clogit(repression_binary ~ issue_performance*v2exl_legitperf + issue_rat_legal*v2exl_legitratio + issue_leader*v2exl_legitlead + issue_regime + topic2 +
            factor(part_violence) + log(avg_participants) + factor(scope) + event_history_location +
            year + strata(cowcode) + cluster(cowcode), data = data, method = "efron")

# get AIC
AIC(m1, m2, m3, m4)

stargazer::stargazer(m1,m2,m3,m4,
                     type = "text",
                     omit = c("cowcode", "year"),
                     digits = 2,
                     star.cutoffs = c(0.05, 0.01, 0.001),
                     covariate.labels = c("Issue: Performance", "Issue: Rational-legal", "Issue: Leader", "Issue: Regime", "Issue: Repression", "Property damage by participants", "People injured by participants", "People killed by participants", "Participants (log)", "Scope: regional/ state", "Scope: local", "Event history", "Claim: Performance", "Claim: Rational-legal", "Claim: Leader", "Performance: Issue x Claim", "Rational-legal: Issue x Claim", "Leader: Issue x Claim"), 
                     dep.var.labels = "",
                     dep.var.caption  = "",
                     column.labels = NULL,
                     #model.numbers = F,
                     #intercept.top = T,
                     #intercept.bottom = F,
                     multicolumn = T,
                     df = F,
                     caption = "Conditional logistic regressions with country- and year-fixed effects. Unit of analysis: All anti-government protest events in MMAD. Robust standard errors clustered at country-level.",           
                     notes.append = T,
                     single.row = F,
                     #object.names = T,
                     style = "apsr",
                     font.size = "footnotesize",
                     add.lines = list(c("Country FEs", "Yes", "Yes", "Yes", "Yes"),
                                      c("Year FEs", "Yes", "Yes", "Yes", "Yes")),
                     title = "")

```

Figures 4a, 4b, and 4c

```{r}
# leader issue and leader claim
d_i1 <- hist(data$v2exl_legitlead, breaks = seq(min(data$v2exl_legitlead),
                                                     max(data$v2exl_legitlead), length.out = 50))
d_x_i1 <- d_i1$breaks[1:49]
d_y_i1 <- d_i1$density

interaction_plot(model = m4, var1 = "issue_leader", var2 = "v2exl_legitlead",
       int = "issue_leader:v2exl_legitlead", vcov = vcov,
       xlab = "Leader legitimation (low-high)", 
       ylab = "Marginal effect of leader issue") + 
  geom_bar(aes(x = d_x_i1, y = d_y_i1), stat = "identity") +
  scale_y_continuous(breaks=c(0,1,2,3,4),
                     labels=c("-2", "-1", "0", "1", "2")) +
    scale_x_continuous(breaks=c(-2,-1,0,1,2,3),
                     labels=c("-2", "-1", "0", "1", "2","3")) +
  coord_cartesian(ylim = c(0, 4)) +
  theme_bw() +
  ggtitle("")


# performance issue and performance claim
d_i2 <- hist(data$v2exl_legitperf, breaks = seq(min(data$v2exl_legitperf),
                                                     max(data$v2exl_legitperf), length.out = 50))
d_x_i2 <- d_i2$breaks[1:49]
d_y_i2 <- d_i2$density

interaction_plot(model = m4, var1 = "issue_performance", var2 = "v2exl_legitperf",
       int = "issue_performance:v2exl_legitperf", vcov = vcov,
       xlab = "Performance legitimation (low-high)", 
       ylab = "Marginal effect of performance issue") + 
  geom_bar(aes(x = d_x_i2, y = d_y_i2), stat = "identity") +
  scale_y_continuous(breaks=c(0,1,2,3,4),
                     labels=c("-2", "-1", "0", "1", "2")) +
  scale_x_continuous(breaks=c(-2,-1,0,1,2,3),
                     labels=c("-2", "-1", "0", "1", "2","3")) +
    coord_cartesian(ylim = c(0, 4)) +
  theme_bw() +
  ggtitle("")


# rational legal issue and rational legal claim
d_i3 <- hist(data$v2exl_legitratio, breaks = seq(min(data$v2exl_legitratio),
                                                     max(data$v2exl_legitratio), length.out = 50))
d_x_i3 <- d_i3$breaks[1:49]
d_y_i3 <- d_i3$density

interaction_plot(model = m4, var1 = "issue_rat_legal", var2 = "v2exl_legitratio",
       int = "issue_rat_legal:v2exl_legitratio", vcov = vcov,
       xlab = "Rational-legal legitimation (low-high)", 
       ylab = "Marginal effect of rational-legal issue") + 
  geom_bar(aes(x = d_x_i3, y = d_y_i3), stat = "identity") +
  scale_y_continuous(breaks=c(0,1,2,3,4),
                     labels=c("-2", "-1", "0", "1", "2")) +
  coord_cartesian(ylim = c(0, 4)) +
  theme_bw() +
  ggtitle("")
```


### Appendix

Table A.1

```{r, echo=FALSE, include=T, eval=T}

# load summary from LDA model
summary_LDA <- read_csv("summary_LDA.csv")

summary_LDA
```

Table A.2

```{r, echo=FALSE, include=T, eval=T}
# calculate pairwise correlations between topisc from BTM and LDA

# repression
cor(data$topic2, data$t_6, use = "pairwise.complete.obs")
# leader
cor(data$topic4, data$t_1, use = "pairwise.complete.obs")
# regime
cor(data$topic7, data$t_2, use = "pairwise.complete.obs")
# elections
cor(data$topic5, data$t_3, use = "pairwise.complete.obs")
# oppositon
cor(data$topic6, data$t_5, use = "pairwise.complete.obs")
# governance
cor(data$topic3, data$t_4, use = "pairwise.complete.obs")
# economy
cor(data$topic8, data$t_7, use = "pairwise.complete.obs")
```


Figure A.1

```{r}
# plot temporal variation of claims for all countries with country-year n > 1

# create helper "mmad" to check for n > 1
mmad <- data %>% select(cowcode, year) %>% mutate(mmad = 1) %>% 
  group_by(cowcode,year) %>% summarise(mmad = mean(mmad))

load("vdem_claims.Rda")

vdem_claims <- vdem_claims  %>% 
  mutate(cowcode = as.factor(cowcode), year = as.factor(year)) %>% 
  left_join(mmad, by = c("cowcode", "year")) %>% 
  filter(!is.na(mmad)) %>% 
  select(-c(mmad,cowcode)) %>% 
  filter(!country_name == "Eritrea") %>% 
  filter(!country_name == "Georgia") %>% 
  filter(!country_name == "Guinea-Bissau") %>% 
  filter(!country_name == "Kenya") %>% 
  filter(!country_name == "Laos") %>% 
  filter(!country_name == "Liberia") %>% 
  filter(!country_name == "Mali") 
  

vdem_claims <- vdem_claims %>% 
  tidyr::pivot_longer(-c(country_name, year)) %>%
  dplyr::mutate(claims = stringr::str_remove(name, "_codelow|_codehigh"),
                type = ifelse(stringr::str_detect(name, "_codelow"), "codelow", "value"),
                type = ifelse(stringr::str_detect(name, "_codehigh"), "codehigh", type)) %>%
  tidyr::pivot_wider(names_from = (type),
                     values_from = value,
                     id_cols = c(country_name, year, claims)) %>%
  tidyr::drop_na() %>% 
  mutate(country_name = if_else(country_name == "Eswatini", "Swaziland", country_name),
         country_name = if_else(country_name == "Democratic Republic of the Congo", "DRC", country_name),
         country_name = if_else(country_name == "Central African Republic", "Central African Rep.", country_name))


ggplot(data = vdem_claims,
       aes(x = year, y = value, group = claims, linetype = claims)) + 
  geom_line()  + 
  geom_ribbon(aes(x = year, ymin = codelow,
                  ymax = codehigh), alpha = .05) +
  xlab("") + ylab("Claim strength") + ylim(0,4) +
  scale_linetype_manual(name = "", values = c(2, 1, 3, 4),
                        labels = c("Ideology", "Leader", "Performance", "Rational-legal")) +
  scale_x_discrete(breaks = seq(2003, 2015, 2)) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white", colour = "black", size = 0.5, linetype = "solid"),
        axis.text.x = element_text(angle = 90, hjust = 1, size = 8),
        panel.grid.major.x = element_line( size=.05, color="grey"),
        legend.position = "bottom",
        strip.text = element_text(size = 8)) +
  facet_wrap(~country_name)
```

Table A.3

```{r}
# summary statistics
summ_stats <- data %>% 
  select("Repression" = repression_binary, "Issue: Performance" = issue_performance, "Issue: Rational-legal" = issue_rat_legal, "Issue: Leader" = issue_leader, "Issue: Regime" = issue_regime, "Issue: Repression"= issue_repression, "Claim: Performance" = v2exl_legitperf, "Claim: Rational-legality" = v2exl_legitratio, "Claim: Leader"= v2exl_legitlead, "Participant violence" = part_violence, avg_participants, "Scope" = scope, "Event history" = event_history_location, "Repressed protests (country, 7 days)" = repr_hist_country_7, "Repressed protests (country, 21 days)" = repr_hist_country_21, "Repressed protests (location, 7 days)" = repr_hist_loc_7, "Repressed protests (location, 21 days)" = repr_hist_loc_21, "Unrepressed protests (country, 7 days)" = unrepr_hist_country_7, "Unrepressed protests (country, 21 days)" = unrepr_hist_country_21, "Unrepressed protests (location, 7 days)" = unrepr_hist_loc_7, "Unrepressed protests (location, 21 days)" = unrepr_hist_loc_21, gdppc, pop, "Electoral democracy index" = v2x_polyarchy, "Respect for human rights" = theta_mean, "Military power base" = v2x_ex_military, "Party power base" = v2x_ex_party, "Armed conflict" = armed_confl, "Personalist power concentration" = latent_personalism) %>% 
  mutate(gdppc = log(gdppc), pop = log(pop), avg_participants = log(avg_participants)) %>% 
  rename("GDP p.c. (log)" = gdppc, "Population (log)" = pop, "Participants (log)" = avg_participants)

summ_stats_table <- psych::describe(summ_stats) %>% select(n, mean, sd, min, max)

xtable(summ_stats_table, digits = 2)
```

Table A.4

```{r}
# with country-level controls
m.a.country <- clogit(repression_binary ~ issue_performance*v2exl_legitperf + issue_rat_legal*v2exl_legitratio + issue_leader*v2exl_legitlead + issue_regime + topic2 +
                        as.factor(part_violence) + log(avg_participants) + as.factor(scope) + event_history_location + 
                        log(gdppc) +  log(pop) + v2x_polyarchy + 
                        theta_mean + v2x_ex_military + v2x_ex_party + armed_confl +
                        year + strata(cowcode) + cluster(cowcode), data = data, method = "efron")

# personalism
m.a.pers <- clogit(repression_binary ~ issue_performance*v2exl_legitperf + issue_rat_legal*v2exl_legitratio + issue_leader*v2exl_legitlead  + issue_regime + topic2 +
                     as.factor(part_violence) + log(avg_participants) + as.factor(scope) + event_history_location + 
                     log(gdppc) +  log(pop) + v2x_polyarchy + theta_mean + v2x_ex_military + v2x_ex_party + armed_confl + latent_personalism + 
                     year + strata(cowcode) + cluster(cowcode), data = data, method = "efron")

# interaction personalism x leader issue
m.a.pers.int <- clogit(repression_binary ~ issue_leader*latent_personalism + issue_performance + v2exl_legitperf + issue_rat_legal + v2exl_legitratio + v2exl_legitlead  + issue_regime + topic2 +
                     as.factor(part_violence) + log(avg_participants) + as.factor(scope) + event_history_location + 
                     log(gdppc) +  log(pop) + v2x_polyarchy + theta_mean + v2x_ex_military + v2x_ex_party + armed_confl +
                     year + strata(cowcode) + cluster(cowcode), data = data, method = "efron")


# table
stargazer::stargazer(m.a.country,m.a.pers,m.a.pers.int,
                     type = "text",
                     omit = c("cowcode", "year"),
                     digits = 2,
                     star.cutoffs = c(0.05, 0.01, 0.001),
                     covariate.labels = c("Issue: Performance", "Claim: Performance", "Issue: Rational-legal", "Claim: Rational-legal", "Issue: Leader", "Claim: Leader", "Issue: Regime", "Issue: Repression", "Property damage by participants", "People injured by participants", "People killed by participants", "Participants (log)", "Scope: regional/ state", "Scope: local", "Event history", "GDP p.c. (log)", "Population (log)", "Electoral democracy index", "Respect for human rights", "Military power base", "Party power base", "Armed conflict", "Personalist power concentration", "Performance: Issue x Claim", "Rational-legal: Issue x Claim", "Leader: Issue x Claim", "Issue: Leader x Personalist power conc."), 
                     #dep.var.labels = c(),
                     dep.var.caption  = "",
                     #model.numbers = F,
                     #intercept.top = T,
                     #intercept.bottom = F,
                     multicolumn = T,
                     df = F,
                     caption = "Conditional logistic regressions. Unit of analysis: Protest events. Robust standard errors clustered at country-level.",           
                     notes.append = T,
                     single.row = F,
                     #object.names = T,
                     style = "apsr",
                     font.size = "footnotesize",
                     add.lines = list(c("Country FEs", "Yes", "Yes", "Yes"),
                                      c("Year FEs", "Yes", "Yes", "Yes")),
                     title = "")
```

Figure A.2

```{r}
# interaction country-level controls
interaction_plot(model = m.a.country, var1 = "issue_leader", var2 = "v2exl_legitlead",
       int = "issue_leader:v2exl_legitlead", vcov = vcov,
       xlab = "Leader legitimation (low-high)", 
       ylab = "Leader issue") + 
  geom_bar(aes(x = d_x_i1, y = d_y_i1), stat = "identity") +
  scale_y_continuous(breaks=c(0,1,2,3,4),
                     labels=c("-2", "-1", "0", "1", "2")) +
  scale_x_continuous(breaks=c(-2,-1,0,1,2,3),
                     labels=c("-2", "-1", "0", "1", "2", "3")) +
  coord_cartesian(ylim = c(0, 4), xlim = c(-2,3)) +
  theme_bw() +
  ggtitle("With country-level controls") +
    theme(plot.title = element_text(size=11))

# interaction country-level controls with personalism
interaction_plot(model = m.a.pers, var1 = "issue_leader", var2 = "v2exl_legitlead",
       int = "issue_leader:v2exl_legitlead", vcov = vcov,
       xlab = "Leader legitimation (low-high)", 
       ylab = "Leader issue") + 
  geom_bar(aes(x = d_x_i1, y = d_y_i1), stat = "identity") +
  scale_y_continuous(breaks=c(0,1,2,3,4),
                     labels=c("-2", "-1", "0", "1", "2")) +
      scale_x_continuous(breaks=c(-2,-1,0,1,2,3),
                     labels=c("-2", "-1", "0", "1", "2", "3")) +
  coord_cartesian(ylim = c(0, 4), xlim = c(-2,3)) +
  theme_bw() +
  ggtitle("With country-level and personalist concentration of power controls") +
  theme(plot.title = element_text(size=11))

```

Figure A.3

```{r}
# interaction personalism
interaction_plot(model = m.a.pers.int, var1 = "issue_leader", var2 = "latent_personalism",
       int = "issue_leader:latent_personalism", vcov = vcov,
       xlab = "Personalist power concentration (low-high)", 
       ylab = "Rational-legal issue") + 
  scale_y_continuous(breaks=c(0,1,2,3,4),
                     labels=c("-2", "-1", "0", "1", "2")) +
  coord_cartesian(ylim = c(0, 4)) +
  theme_bw() +
  ggtitle("")
```

Table A.5

```{r}
# claims
m.a.maxclaim.perf <- clogit(repression_binary ~  issue_performance*max_claim_perf + issue_regime + issue_repression +
                        as.factor(part_violence) + log(avg_participants) + as.factor(scope) + event_history_location + 
                        log(gdppc) +  log(pop) + v2x_polyarchy + theta_mean + v2x_ex_military + v2x_ex_party + armed_confl +
                        year + strata(cowcode) + cluster(cowcode), data = data, method = "efron")

m.a.maxclaim.ratio <- clogit(repression_binary ~  issue_rat_legal*max_claim_ratio + issue_regime + issue_repression +
                        as.factor(part_violence) + log(avg_participants) + as.factor(scope) + event_history_location + 
                        log(gdppc) +  log(pop) + v2x_polyarchy + theta_mean + v2x_ex_military + v2x_ex_party + armed_confl +
                        year + strata(cowcode) + cluster(cowcode), data = data, method = "efron")

m.a.maxclaim.lead <- clogit(repression_binary ~  issue_leader*max_claim_lead + issue_regime + issue_repression +
                        as.factor(part_violence) + log(avg_participants) + as.factor(scope) + event_history_location + 
                        log(gdppc) +  log(pop) + v2x_polyarchy + theta_mean + v2x_ex_military + v2x_ex_party + armed_confl +
                        year + strata(cowcode) + cluster(cowcode), data = data, method = "efron")

# ideology as reference category
m.a.maxclaim.all <- clogit(repression_binary ~  issue_performance*max_claim_perf + issue_rat_legal*max_claim_ratio + issue_leader*max_claim_lead + issue_regime + issue_repression +
                        as.factor(part_violence) + log(avg_participants) + as.factor(scope) + event_history_location + 
                        log(gdppc) +  log(pop) + v2x_polyarchy + theta_mean + v2x_ex_military + v2x_ex_party + armed_confl +
                        year + strata(cowcode) + cluster(cowcode), data = data, method = "efron")


# table
stargazer::stargazer(m.a.maxclaim.perf, m.a.maxclaim.ratio, m.a.maxclaim.lead, m.a.maxclaim.all,
                     type = "text", 
                     omit = c("cowcode","year"),
                     digits = 2,
                     star.cutoffs = c(0.05, 0.01, 0.001),
                     covariate.labels = c("Issue: Performance","Claim (max.): Performance", "Issue: Rational-legal", "Claim (max.): Rational-legal",  "Issue: Leader",  "Claim (max.): Leader", "Issue: Regime", "Issue: Repression", "Property damage by participants", "People injured by participants", "People killed by participants","Participants (log)", "Scope: regional/ state", "Scope: local", "Event history","GDP p.c. (log)", "Population (log)", "Electoral democracy index", "Respect for human rights", "Military power base", "Party power base", "Armed conflict", "Performance: Issue x Claim (max.)", "Rational-legal: Issue x Claim (max.)", "Leader: Issue x Claim (max.)"), 
                     #dep.var.labels = c(),
                     dep.var.caption  = "",
                     #model.numbers = F,
                     #intercept.top = T,
                     #intercept.bottom = F,
                     multicolumn = T,
                     df = F,
                     caption = "Conditional logistic regressions. Unit of analysis: Protest events. Robust standard errors clustered at country-level. 'Claim: Ideology (max.)' as reference category in Model 11.",
                     notes.append = T,
                     single.row = F,
                     #object.names = T,
                     style = "apsr",
                     font.size = "footnotesize",
                     add.lines = list(c("Country FEs", "Yes", "Yes", "Yes", "Yes"),
                                      c("Year FEs", "Yes", "Yes", "Yes", "Yes")),
                     label = "tab:tab_app_max_claim")

```

Table A.6

```{r}
# issues
m.a.maxissue.perf <- clogit(repression_binary ~ max_issue_perf*v2exl_legitperf + max_issue_regime + 
                        as.factor(part_violence) + log(avg_participants) + as.factor(scope) + event_history_location + 
                        log(gdppc) +  log(pop) + v2x_polyarchy + theta_mean + v2x_ex_military + v2x_ex_party + armed_confl +
                        year + strata(cowcode) + cluster(cowcode), data = data, method = "efron")

m.a.maxissue.ratio <- clogit(repression_binary ~ max_issue_ratio*v2exl_legitratio + max_issue_regime +
                        as.factor(part_violence) + log(avg_participants) + as.factor(scope) + event_history_location + 
                        log(gdppc) +  log(pop) + v2x_polyarchy + theta_mean + v2x_ex_military + v2x_ex_party + armed_confl +
                        year + strata(cowcode) + cluster(cowcode), data = data, method = "efron")

m.a.maxissuelead <- clogit(repression_binary ~ max_issue_lead*v2exl_legitlead  + max_issue_regime +
                        as.factor(part_violence) + log(avg_participants) + as.factor(scope) + event_history_location + 
                        log(gdppc) +  log(pop) + v2x_polyarchy + theta_mean + v2x_ex_military + v2x_ex_party + armed_confl +
                        year + strata(cowcode) + cluster(cowcode), data = data, method = "efron")

m.a.maxissue.all <- clogit(repression_binary ~ max_issue_perf*v2exl_legitperf + max_issue_ratio*v2exl_legitratio + max_issue_lead*v2exl_legitlead  + max_issue_regime +
                        as.factor(part_violence) + log(avg_participants) + as.factor(scope) + event_history_location + 
                        log(gdppc) +  log(pop) + v2x_polyarchy + theta_mean + v2x_ex_military + v2x_ex_party + armed_confl +
                        year + strata(cowcode) + cluster(cowcode), data = data, method = "efron")

stargazer::stargazer(m.a.maxissue.perf, m.a.maxissue.ratio, m.a.maxissuelead, m.a.maxissue.all, 
                     type = "text", 
                     omit = c("cowcode","year"),
                     digits = 2,
                     star.cutoffs = c(0.05, 0.01, 0.001),
                     covariate.labels = c("Issue (max.): Performance","Claim: Performance", "Issue (max.): Rational-legal", "Claim: Rational-legal",  "Issue (max): Leader",  "Claim: Leader", "Issue (max.): Regime" , "Property damage by participants", "People injured by participants", "People killed by participants","Participants (log)", "Scope: regional/ state", "Scope: local", "Event history","GDP p.c. (log)", "Population (log)", "Electoral democracy index", "Respect for human rights", "Military power base", "Party power base", "Armed conflict", "Performance: Issue (max.) x Claim", "Rational-legal: Issue (max.) x Claim", "Leader: Issue (max.) x Claim"), 
                     #dep.var.labels = c(),
                     dep.var.caption  = "",
                     #model.numbers = F,
                     #intercept.top = T,
                     #intercept.bottom = F,
                     multicolumn = T,
                     df = F,
                     caption = "Conditional logistic regressions. Unit of analysis: Protest events. Robust standard errors clustered at country-level.", 
                     notes.append = T,
                     single.row = F,
                     #object.names = T,
                     style = "apsr",
                     font.size = "footnotesize",
                     add.lines = list(c("Country FEs", "Yes", "Yes", "Yes", "Yes"),
                                      c("Year FEs", "Yes", "Yes", "Yes", "Yes")),
                     label = "tab:tab_app_max_issue")

```

Figure A.4

```{r}
interaction_plot_binary(model = m.a.maxclaim.all, var1 = "issue_leader", var2 = "max_claim_lead",
                                     int = "issue_leader:max_claim_lead", vcov = vcov,
                                     xlab = "Leader legitimation strongest claim", 
                                     ylab = "Marginal effect of leader issue") +
  scale_x_continuous(breaks = c(0,1), labels = c("no", "yes")) +
  scale_y_continuous(breaks=c(-0.5,-0.25,0.0,0.25,0.5,0.75,1.00,1.25), labels=c("-0.50","-0.25", "0.00", "0.25", "0.50", "0.75", "1.00", "1.25")) +
  theme_bw() +
  ggtitle("")
```

Figure A.5

```{r}
interaction_plot(model = m.a.maxissue.all, var1 = "max_issue_lead", var2 = "v2exl_legitlead",
                 int = "max_issue_lead:v2exl_legitlead", vcov = vcov,
                 xlab = "Leader legitimation (low-high)", 
                 ylab = "Leader strongest issue") + 
  geom_bar(aes(x = d_x_i1, y = d_y_i1), stat = "identity") +
scale_y_continuous(breaks=c(0,1,2,3,4,5),
                    labels=c("-2", "-1", "0", "1", "2","3")) +
  coord_cartesian(ylim = c(0, 4), xlim = c(-2,3)) +
  theme_bw() +
  ggtitle("")

```

Table A.7

```{r}
m.a.repr.hist.c.7 <- clogit(repression_binary ~ issue_performance*v2exl_legitperf + repr_hist_country_7 + unrepr_hist_country_7 + issue_rat_legal*v2exl_legitratio + issue_leader*v2exl_legitlead + issue_regime + topic2 +
            factor(part_violence) + log(avg_participants) + factor(scope) +
            year + strata(cowcode) + cluster(cowcode), data = data, method = "efron")

m.a.repr.hist.c.21 <- clogit(repression_binary ~ issue_performance*v2exl_legitperf + repr_hist_country_21 + unrepr_hist_country_21 + issue_rat_legal*v2exl_legitratio + issue_leader*v2exl_legitlead + issue_regime + topic2 +
            factor(part_violence) + log(avg_participants) + factor(scope) + 
            year + strata(cowcode) + cluster(cowcode), data = data, method = "efron")

m.a.repr.hist.l.7 <- clogit(repression_binary ~ issue_performance*v2exl_legitperf + repr_hist_loc_7 + unrepr_hist_loc_7 + issue_rat_legal*v2exl_legitratio + issue_leader*v2exl_legitlead + issue_regime + topic2 +
            factor(part_violence) + log(avg_participants) + factor(scope) +
            year + strata(cowcode) + cluster(cowcode), data = data, method = "efron")

m.a.repr.hist.l.21 <- clogit(repression_binary ~ issue_performance*v2exl_legitperf + repr_hist_loc_21 + unrepr_hist_loc_21 + issue_rat_legal*v2exl_legitratio + issue_leader*v2exl_legitlead + issue_regime + topic2 +
            factor(part_violence) + log(avg_participants) + factor(scope) + 
            year + strata(cowcode) + cluster(cowcode), data = data, method = "efron")

stargazer::stargazer(m.a.repr.hist.c.7, m.a.repr.hist.c.21, m.a.repr.hist.l.7, m.a.repr.hist.l.21,
                     type = "text", 
                     omit = c("cowcode", "year"), 
                     star.cutoffs = c(0.05, 0.01, 0.001),
                     digits = 2,
                     covariate.labels = c("Issue: Performance", "Claim: Performance", "Repressed protests (country, 7 days)", "Unrepressed protest (country, 7 days)", "Repressed protests (country, 21 days)", "Unrepressed protest (country, 21 days)", "Repressed protests (location, 7 days)", "Unrepressed protest (location, 7 days)", "Repressed protests (location, 21 days)", "Unrepressed protest (location, 21 days)", "Issue: Rational-legal", "Claim: Rational-legal", "Issue: Leader",  "Claim: Leader", "Issue: Regime", "Issue: Repression", "Property damage by participants", "People injured by participants", "People killed by participants", "Participants (log)", "Scope: regional/ state", "Scope: local", "Performance: Issue x Claim", "Rational-legal: Issue x Claim", "Leader: Issue x Claim"), 
                     dep.var.labels = "",
                     dep.var.caption  = "",
                     column.labels = NULL,
                     multicolumn = T,
                     df = F,
                     caption = "Conditional logistic regressions with country- and year-fixed effects. Unit of analysis: All anti-government protest events in MMAD. Robust standard errors clustered at country-level.",           
                     notes.append = T,
                     single.row = F,
                     style = "apsr",
                     font.size = "footnotesize",
                     add.lines = list(c("Country FEs", "Yes", "Yes", "Yes", "Yes"),
                                      c("Year FEs", "Yes", "Yes", "Yes", "Yes")),
                     title = "",
                     label = "tab:tab_a_repr_history")
```

Figure A.6

```{r}
interaction_plot(model = m.a.repr.hist.c.7, var1 = "issue_leader", var2 = "v2exl_legitlead",
       int = "issue_leader:v2exl_legitlead", vcov = vcov,
       xlab = "Leader legitimation (low-high)", 
       ylab = "Marginal effect of leader issue") + 
  geom_bar(aes(x = d_x_i1, y = d_y_i1), stat = "identity") +
  scale_y_continuous(breaks=c(0,1,2,3,4),
                     labels=c("-2", "-1", "0", "1", "2")) +
    scale_x_continuous(breaks=c(-2,-1,0,1,2,3),
                     labels=c("-2", "-1", "0", "1", "2","3")) +
  coord_cartesian(ylim = c(0, 4)) +
  theme_bw() +
  ggtitle("")
```

Figure A.7

```{r}
interaction_plot(model = m.a.repr.hist.l.7, var1 = "issue_leader", var2 = "v2exl_legitlead",
       int = "issue_leader:v2exl_legitlead", vcov = vcov,
       xlab = "Leader legitimation (low-high)", 
       ylab = "Marginal effect of leader issue") + 
  geom_bar(aes(x = d_x_i1, y = d_y_i1), stat = "identity") +
  scale_y_continuous(breaks=c(0,1,2,3,4),
                     labels=c("-2", "-1", "0", "1", "2")) +
    scale_x_continuous(breaks=c(-2,-1,0,1,2,3),
                     labels=c("-2", "-1", "0", "1", "2","3")) +
  coord_cartesian(ylim = c(0, 4)) +
  theme_bw() +
  ggtitle("")
```

Table A.9

```{r}
# aggregate events to yearly obs
df <- data %>%
    dplyr::group_by(cowcode, year) %>% 
    dplyr::summarise(events = n(), gdppc = mean(gdppc), pop = mean(pop), v2exl_legitlead = mean(v2exl_legitlead), v2exl_legitperf = mean(v2exl_legitperf), v2exl_legitratio = mean(v2exl_legitratio), v2x_polyarchy = mean(v2x_polyarchy), theta_mean = mean(theta_mean), issue_performance = mean(issue_performance), issue_leader = mean(issue_leader), issue_rat_legal = mean(issue_rat_legal), issue_regime = mean(issue_regime)) %>% 
    ungroup()

# yearly number of events
m.a.events <- glm(events ~ v2exl_legitlead + v2exl_legitperf + v2exl_legitratio + log(gdppc) + log(pop)+ v2x_polyarchy + theta_mean + as.factor(year) + as.factor(cowcode), family = "poisson", data = df)

m.a.events.se <- coeftest(m.a.events, cluster.vcov(m.a.events, df$cowcode))

# issue strength
m.a.iss.perf <- lm(issue_performance ~ v2exl_legitperf + v2exl_legitlead + v2exl_legitratio + log(gdppc) + log(pop) + v2x_polyarchy + theta_mean + cowcode + year, data = data)

m.a.iss.perf.se <- coeftest(m.a.iss.perf, vcov = vcovHC(m.a.iss.perf,type="HC0"), cluster = ~ cowcode)

m.a.iss.lead <- lm(issue_leader ~ v2exl_legitperf + v2exl_legitlead + v2exl_legitratio + log(gdppc) + log(pop) + v2x_polyarchy + theta_mean + cowcode + year, data = data)

m.a.iss.lead.se <- coeftest(m.a.iss.lead, vcov = vcovHC(m.a.iss.lead,type="HC0"), cluster = ~ cowcode)

m.a.iss.rat <- lm(issue_rat_legal~ v2exl_legitperf + v2exl_legitlead + v2exl_legitratio + log(gdppc) + log(pop) + v2x_polyarchy + theta_mean + cowcode + year, data = data)

m.a.iss.rat.se <- coeftest(m.a.iss.rat, vcov = vcovHC(m.a.iss.rat,type="HC0"), cluster = ~ cowcode)

stargazer::stargazer(m.a.events.se, m.a.iss.perf.se, m.a.iss.lead.se, m.a.iss.rat.se,
    type = "text",
    omit = c("cowcode", "year"),
    digits = 2,
    star.cutoffs = c(0.05, 0.01, 0.001),
    covariate.labels = c("Claim: Leader", "Claim: Performance", "Claim: Rational-legal", "GDP p.c. (log)", "Population (log)", "Electoral democracy index", "Respect for human rights"),
    column.labels = c("Yearly number of events", "Issue: Performance", "Issue: Leader", "Issue: Rational-legal"),
    dep.var.caption  = "",
    multicolumn = T,
          df = F,
          caption = "Model 22: Poisson regressions. DV: Yearly number of protest events. Unit of observation: Country-year. Models 23-25: OLS regressions. DVs: Issue strengths. Unit of analysis: All anti-government protest events in MMAD. Robust standard errors clustered by country.",
          notes.append = T,
          single.row = F,
          #object.names = T,
          style = "apsr",
          font.size = "footnotesize",
          add.lines = list(c("Country FEs", "Yes", "Yes", "Yes", "Yes"),
                           c("Year FEs","Yes", "Yes", "Yes", "Yes"),
                           c("N", "601", "14,742", "14,742", "14,742")),
          label = "tab:tab_app_events_issues")
```

Figure A.9

```{r}
# select strongest claim, get event number per country-year
data.agg <- data %>% 
  dplyr::group_by(cowcode,year, max_claim) %>% 
  summarise(events = length(event_id)) %>% 
  ungroup() 

data.agg$max_claim <- as.factor(data.agg$max_claim)

# plot number of events per year and country across strongest claims
ggplot(data.agg, aes(events, group = max_claim)) +
  geom_density(aes(linetype = max_claim), show_guide = FALSE)+
  stat_density(aes(x = events, linetype = max_claim), 
               geom = "line", position = "identity") +
  theme_classic() +
  theme(legend.title = element_blank()) +
  labs(title = NULL,
       x = NULL,
       y = NULL,
       fill = "") +
  xlim(0,100) +
  theme(legend.position="bottom", legend.spacing.x = unit(0.1, 'cm')) +
  scale_linetype_manual(labels = c("Ideology", "Leader", "Performance", "Rational-legality"),
                        values = c("dotdash", "dashed", "dotted", "solid"))

```

Figure A.10

```{r}
# plot issues across claims, selecting events with above average strength
mean(data$issue_performance, na.rm = T)

data %>% 
  filter(issue_performance > 0.1929606) %>% 
  pivot_longer(cols = c(v2exl_legitlead, v2exl_legitperf, v2exl_legitratio), names_to = "claim",
               values_to = "strength") %>% 
  ggplot() +
  geom_density(aes(x = strength, linetype = claim), show_guide = FALSE)+
  stat_density(aes(x = strength, linetype = claim), 
               geom = "line", position = "identity") +
  theme_classic() +
  theme(legend.title = element_blank(),
        plot.title = element_text(size = 13),
        legend.spacing.x = unit(0.1, 'cm'),
        legend.text=element_text(size = 13),
        axis.text=element_text(size = 13)) +
  labs(title = "Protest raising above average performance issues",
       x = NULL,
       y = NULL,
       fill = "") +
  xlim(-3,3) +
  ylim(0,0.8) +
  scale_linetype_manual(labels = c("Leader legitimation", "Performance legitimation",
                                   "Rational-legal legitimation"),
                        values = c("dashed", "dotted", "solid"))

mean(data$issue_rat_legal, na.rm = T)

data %>% 
  filter(issue_rat_legal > 0.2309849) %>% 
  pivot_longer(cols = c(v2exl_legitlead, v2exl_legitperf, v2exl_legitratio), names_to = "claim",
               values_to = "strength") %>% 
  ggplot() +
  geom_density(aes(x = strength, linetype = claim), show_guide = FALSE)+
  stat_density(aes(x = strength, linetype = claim), 
               geom = "line", position = "identity") +
  theme_classic() +
  theme(legend.title = element_blank(),
        plot.title = element_text(size = 13),
        legend.spacing.x = unit(0.1, 'cm'),
        axis.text=element_text(size = 13)) +
  labs(title = "Protest raising above average rational-legal issues",
       x = NULL,
       y = NULL,
       fill = "") +
  theme(legend.position="none") +
  xlim(-3,3) +
  ylim(0,0.8) 

mean(data$issue_leader, na.rm = T)

data %>% 
  filter(issue_leader > 0.1235358) %>%  
  pivot_longer(cols = c(v2exl_legitlead, v2exl_legitperf, v2exl_legitratio), names_to = "claim",
               values_to = "strength") %>% 
  ggplot() +
  geom_density(aes(x = strength, linetype = claim), show_guide = FALSE)+
  stat_density(aes(x = strength, linetype = claim), 
               geom = "line", position = "identity") +
  theme_classic() +
  theme(legend.title = element_blank(),
        legend.position = "bottom",
        plot.title = element_text(size = 13),
        legend.spacing.x = unit(0.1, 'cm'),
        axis.text=element_text(size = 13)) +
  labs(title = "Protest raising above average leader issues",
       x = NULL,
       y = NULL,
       fill = "") +
  theme(legend.position="none") +
  xlim(-3,3) +
  ylim(0,0.8) 

```
